﻿// h8truecolor_yanjiuyuan.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
// 给研究院做的葵花8 h8 真彩图程序 2020-8-31


#include "gdal_priv.h"

#include <string>
#include <iostream>
#include <vector>
#include <Windows.h>
#include <cmath>
//#include "../sharedcodes2/wftools.h"
//#include "../sharedcodes2/ajson5.h"
#include <cassert>
#include "H5Cpp.h"

#include <cstdio>
#include <fstream>
 
#include <algorithm>
#define WGDALRASTER_H_IMPLEMENTATION
#include "wGdalRaster.h"
 
 


//using namespace H5;
using namespace std;
//using namespace ArduinoJson;

#define DEG2RAD	0.0174532925199		/* PI/180 */
#define UO3	0.319
#define UH2O	2.93
#define REFLMIN -0.01
#define REFLMAX  1.6

#define MAXSOLZ 86.5
#define MAXAIRMASS 28
#define	SCALEHEIGHT 8000
#define MAXNUMSPHALBVALUES	3000
#define	TAUMAX	0.3
#define	TAUSTEP4SPHALB	(TAUMAX / (float)MAXNUMSPHALBVALUES)
 

void sunpos(
	int iYear,
	int iMonth,
	int iDay,
	int dHours,
	int dMinutes,
	int dSeconds,
	double dLongitude,
	double dLatitude,
	double& dZenithAngle,
	double& dAzimuth);


double fintexp1(float tau)
{
	double xx, xftau;
	int i;
	const float a[6] = { -.57721566, 0.99999193,-0.24991055,
					0.05519968,-0.00976004, 0.00107857 };
	xx = a[0];
	xftau = 1.;
	for (i = 1; i < 6; i++) {
		xftau *= tau;
		xx += a[i] * xftau;
	}
	return xx - logf(tau);
}

//! to compute exponential integrads for argument tau. see wiki for detail.
double fintexp3(float tau)
{
	return (expf(-tau) * (1. - tau) + tau * tau * fintexp1(tau)) / 2.;
}

//! to compute the spherical albedo of molecular layer. see 6s manual part2 page96.
/// 计算大气球反照率
/// @param tau Rayleigh光学厚度
/// @return 返回大气球反照率
float csalbr(float tau)
{
	return (3 * tau - fintexp3(tau) * (4 + 2 * tau) + 2 * expf(-tau)) / (4 + 3 * tau);
}


void chand(
	float phi,
	float muv,
	float mus,
	float* taur,
	float* rhoray,
	float* trup,
	float* trdown
)
{
	const double xfd = 0.958725775;
	const float xbeta2 = 0.5;
	float pl[5];
	double fs01, fs02, fs0, fs1, fs2;
	const float as0[10] = { 0.33243832, 0.16285370, -0.30924818, -0.10324388, 0.11493334,
						   -6.777104e-02, 1.577425e-03, -1.240906e-02, 3.241678e-02, -3.503695e-02 };
	const float as1[2] = { .19666292, -5.439061e-02 };
	const float as2[2] = { .14545937,-2.910845e-02 };
	float phios, xcosf1, xcosf2, xcosf3;
	float xph1, xph2, xph3, xitm1, xitm2;
	float xlntaur, xitot1, xitot2, xitot3;
	int i, ib;

	phios = phi + 180;
	xcosf1 = 1.;
	xcosf2 = cosf(phios * DEG2RAD);
	xcosf3 = cosf(2 * phios * DEG2RAD);
	xph1 = 1 + (3 * mus * mus - 1) * (3 * muv * muv - 1) * xfd / 8.;
	xph2 = -xfd * xbeta2 * 1.5 * mus * muv * sqrtf(1 - mus * mus) * sqrtf(1 - muv * muv);
	xph3 = xfd * xbeta2 * 0.375 * (1 - mus * mus) * (1 - muv * muv);
	pl[0] = 1.;
	pl[1] = mus + muv;
	pl[2] = mus * muv;
	pl[3] = mus * mus + muv * muv;
	pl[4] = mus * mus * muv * muv;
	fs01 = fs02 = 0;
	for (i = 0; i < 5; i++) fs01 += pl[i] * as0[i];
	for (i = 0; i < 5; i++) fs02 += pl[i] * as0[5 + i];
	for (ib = 0; ib < 4; ib++) {
		xlntaur = logf(taur[ib]);
		fs0 = fs01 + fs02 * xlntaur;
		fs1 = as1[0] + xlntaur * as1[1];
		fs2 = as2[0] + xlntaur * as2[1];
		trdown[ib] = expf(-taur[ib] / mus);
		trup[ib] = expf(-taur[ib] / muv);
		xitm1 = (1 - trdown[ib] * trup[ib]) / 4. / (mus + muv);
		xitm2 = (1 - trdown[ib]) * (1 - trup[ib]);
		xitot1 = xph1 * (xitm1 + xitm2 * fs0);
		xitot2 = xph2 * (xitm1 + xitm2 * fs1);
		xitot3 = xph3 * (xitm1 + xitm2 * fs2);
		rhoray[ib] = xitot1 * xcosf1 + xitot2 * xcosf2 * 2 + xitot3 * xcosf3 * 2;
	}
}


int getatmvariables(
	float mus,
	float muv,
	float phi,
	short height,
	float* sphalb,
	float* rhoray,
	float* TtotraytH2O,
	float* tOG
)
{
	double m, Ttotrayu, Ttotrayd, tO3, tO2, tH2O, psurfratio;
	int j, ib;

	
	///
	const float aH2O[4] = { 0 ,           0  ,   -5.25251 , -5.25251 };//b g r nir
	const float bH2O[4] = { 0 ,           0  ,   0.820175 , 0.725159 };// b g r  nir 
	const float aO3[4] = {  0.007432 , 0.08969 , 0.07153  , 0.000 };//b g r nir 

	// 6s for H8 // const float taur0[4] = { 0.1854 , 0.135  , 0.0542  ,  0.01629 }; //blue green red nir
	const float taur0[4] = { 0.1630 , 0.112  , 0.04600  ,  0.01629 };

	float taur[4], trup[4], trdown[4];
	static float sphalb0[MAXNUMSPHALBVALUES];
	static bool first_time = true;

	if (first_time) {
		sphalb0[0] = 0;
		for (j = 1; j < MAXNUMSPHALBVALUES; j++)		/* taur <= 0.3 for bands 1 to 7 (including safety margin for height<~0) */
			sphalb0[j] = csalbr(j * TAUSTEP4SPHALB);//! compute molcular spherical albedo from tau 0.0 to 0.3.
		first_time = false;
	}

	m = 1 / mus + 1 / muv;
	if (m > MAXAIRMASS) return -1;
	psurfratio = expf(-height / (float)SCALEHEIGHT);
	for (ib = 0; ib <4 ; ib++)
		taur[ib] = taur0[ib] * psurfratio;

	chand(phi, muv, mus, taur, rhoray, trup, trdown);

	for (ib = 0; ib < 4; ib++) {
		int tauindex = (int)(taur[ib] / TAUSTEP4SPHALB + 0.5); //bug fixed for tau greater 0.3
		tauindex = min(MAXNUMSPHALBVALUES - 1, tauindex);
		sphalb[ib] = sphalb0[tauindex];
		Ttotrayu = ((2 / 3. + muv) + (2 / 3. - muv) * trup[ib]) / (4 / 3. + taur[ib]);
		Ttotrayd = ((2 / 3. + mus) + (2 / 3. - mus) * trdown[ib]) / (4 / 3. + taur[ib]);
		tO3 = tO2 = tH2O = 1;
		if (aO3[ib] != 0) tO3 = expf(-m * UO3 * aO3[ib]);
		if (bH2O[ib] != 0) tH2O = expf(-expf(aH2O[ib] + bH2O[ib] * logf(m * UH2O)));
		TtotraytH2O[ib] = Ttotrayu * Ttotrayd * tH2O;
		tOG[ib] = tO3 * tO2;
	}
	return 0;
}




/// 针对特定波段进行快速大气校正计算
/// @param refl 该波段卫星接收到的大气顶层反射率[0~1.0]
/// @param TtotraytH2O 该波段Rayleigh和H2O综合透过率
/// @param tOG 该波段O2和O3综合透过率
/// @param rhoray 该波段Rayleigh散射反射率
/// @param sphalb 该波段球反照率
/// @return 返回大气校正后的反射率
float correctedrefl(float refl, float TtotraytH2O, float tOG, float rhoray, float sphalb)
{
	float corr_refl;
	//corr_refl = (refl / tOG - rhoray) / TtotraytH2O;
	corr_refl = (refl / tOG - rhoray) / TtotraytH2O;
	corr_refl /= (1 + corr_refl * sphalb);
	return corr_refl;
}

int linear255(double val)
{
	if (val < 0) val = 0;
	return MIN(val * 231.9f, 255);
}

int jac255(int val)
{
	if (val < 30)
	{
		return  val * 3.667f;
	}
	else if (val < 60)
	{
		return val * 1.667f + 60;
	}
	else if (val < 120)
	{
		return val * 0.833f + 110;
	}
	else if (val < 190)
	{
		return val * 0.4286f + 158.57;
	}
	else
	{
		return MIN(255, val * 0.2308f + 196.15);
	}
}

string getDirByFilepath(string str)
{
	size_t found = str.find_last_of("/\\");
	if (found == string::npos) return "";
	else
	{
		return str.substr(0, found + 1);
	}
}

/// 初始化GDAL路径配置相关环境变量
///2018-3-16
void initGdal(string& exeDir)
{
	//cout<<"gdal-driver-plugin:"<<putenv("GDAL_DRIVER_PATH=E:/coding/fy4qhzx-project/fy4prepare/x64/Release")<<endl ;
	//cout<<"proj:"<<putenv("PROJ_LIB=\\proj")<<endl ;

	/// @note 程序必须在运行前设置运行目录，如果程序已经启动那么在设置运行路径也无法找到dll
	/// @note 系统调用的套路是exe启动的同时检查dll是否完整，因此进入程序后无法再设置工作目录，所以dll只能事前放在exe目录下。

	char  szPath[1024];
	GetModuleFileNameA(NULL, szPath, 1024);
	cout << "module:" << szPath << endl;
	string dir = getDirByFilepath(string(szPath)) + "";
	cout << "dir:" << dir << endl;
	exeDir = dir;

	string gdaldataStr = string("GDAL_DATA=") + dir + "data";
	string dataDirStr = dir + "data";
	CPLSetConfigOption("GDAL_DATA", dataDirStr.c_str());
	// cout<<"gdal-data:"<<gdaldataStr<<" putenv:"<<putenv(gdaldataStr.c_str())<<endl ;

	string pluginStr = string("GDAL_DRIVER_PATH=") + dir;
	string driverDir = dir;
	CPLSetConfigOption("GDAL_DRIVER_PATH", driverDir.c_str());
	// cout<<"gdal-driver-plugin:"<<pluginStr<<" putenv:"<<putenv(pluginStr.c_str())<<endl ;

	//string projStr = string("PROJ_LIB=") + dir + "proj" ;
	//cout<<"proj:"<<projStr<<" putenv:"<<putenv(projStr.c_str())<<endl ;

	//cout<<"change working dir:"<<dir<<endl ;
	//string dirstr = "chdir " + dir ;
	//_chdir( dirstr.c_str() ) ;
	// system( dirstr.c_str() ) ; //这个命令会弹出命令行窗口
}
int main(int argc, char* argv[])
{

	cout << "A testing program to make h8 geocolor image. 2020-8-31" << endl;
	cout << "call:h8truecolor_nanjin inputh8.hdf leftlon toplat resolution dem.tif output.tif" << endl;
	cout << "for 6001x6001 data: left 80.0 top 60.0 reso 0.02. " << endl;
	cout << "v1.0" << endl;
	cout << "v1.0.1 upload to gitee, add extra info." << endl;
	cout << "v1.0.2 add geotrans and proj." << endl;
	string exedir;
	initGdal(exedir);
	if (argc != 7 )
	{
		cout << "params number not equal 7." << endl;
		exit(11);
	}

	GDALAllRegister();

	//load h8
	string h8file = argv[1];
	string blueband = "HDF5:" + h8file + "://albedo_01";//针对王金英数据适配
	string greenband = "HDF5:" + h8file + "://albedo_02";
	string redband = "HDF5:" + h8file + "://albedo_03";
	string nirband = "HDF5:" + h8file + "://albedo_04";
	string u10band = "HDF5:" + h8file + "://tbb_13";//13
	string u03band = "HDF5:" + h8file + "://tbb_07";//7

	string t0file = "HDF5:" + h8file + "://SOZ";
	string p0file = "HDF5:" + h8file + "://SOA";
	string t1file = "HDF5:" + h8file + "://SAZ";
	string p1file = "HDF5:" + h8file + "://SAA";

	//left top longlat
	double leftLon = atof(argv[2]);
	double topLat = atof(argv[3]);
	double resolution = atof(argv[4]);
	

	//load dem
	string demfile = argv[5];


	//output file
	string outfile = argv[6];


	wGdalRasterU16 bluedata, greendata, reddata,nirdata , u10data,u03data , t0data, t1data,p0data,p1data ;
	wGdalRasterI16 demdata;
	bool isokblue = bluedata.open(blueband);
	if (isokblue == false)
	{
		cout << "Error : failed to open blue " << blueband << endl;
		return 13;
	}
	greendata.open(greenband);
	reddata.open(redband);
	nirdata.open(nirband);
	u10data.open(u10band);
	u03data.open(u03band);

	t0data.open(t0file);
	t1data.open(t1file);
	p0data.open(p0file);
	p1data.open(p1file);

	bool demok = demdata.open(demfile);
	if (demok == false)
	{
		cout << "Error : failed to open " << demfile << endl;
		return 14;
	}
	const double* demTrans = demdata.getTrans();

	size_t xsize = bluedata.getXSize();
	size_t ysize = bluedata.getYSize();

	wGdalRasterByte outdata;
	outdata.create(xsize, ysize, 3);

	int per0 = -1;

	double gamma = 2;
	double gammainv = 1 / gamma;
	for (size_t iy = 0; iy < ysize; ++iy)
	{
		for (size_t ix = 0; ix < xsize; ++ix)
		{
			int bdn = bluedata.getValuei(ix, iy, 0);
			int gdn = greendata.getValuei(ix, iy, 0);
			int rdn = reddata.getValuei(ix, iy, 0);
			int ndn = nirdata.getValuei(ix, iy, 0);
			int u10dn = u10data.getValuei(ix, iy, 0);
			int u03dn = u03data.getValuei(ix, iy, 0);
			int t0 = t0data.getValuei(ix, iy, 0);
			int t1 = t1data.getValuei(ix, iy, 0);
			int p0 = p0data.getValuei(ix, iy, 0);
			int p1 = p1data.getValuei(ix, iy, 0);

			if (bdn >= 0 && bdn < 15000 &&   u10dn < 60000  && t1 < 8900 )
			{//good pixel
				//night color
				int nightr = 0;
				int nightg = 0;
				int nightb = 0;

				//ice cloud
				int u10u3diff = u10dn - u03dn;
				double u10tbb = u10dn / 100.0 + 273.15;//band13

				if (u10u3diff > 30)
				{
					nightr = MIN(255,MAX(0,0.363 * u10u3diff + 1.215));
					nightg = MIN(255, MAX(0, 0.490 * u10u3diff + 0.1113));
					nightb = MIN(255, MAX(0, 0.640 * u10u3diff - 0.5526));
				}
				else
				{
					//ir cloud
					
					int nightGray = u10tbb * (-2.451) + 688;
					if (nightGray < 0) nightGray = 0;
					else if (nightGray > 255) nightGray = 255;
					nightr = nightGray;
					nightg = nightGray;
					nightb = nightGray;
				}

				if (t0 < 8500)
				{//daytime
					
					double bref = bdn / 10000.0;
					double gref = gdn / 10000.0;
					double rref = rdn / 10000.0;
					double nref = ndn / 10000.0;
					
					//angles
					double t0d = t0 / 100.0;
					double t1d = t1 / 100.0;
					double p0d = p0 / 100.0;
					double p1d = p1 / 100.0;
					double mu0 = cos(t0d * DEG2RAD);
					double mu1 = cos(t1d * DEG2RAD);
					double rph = abs(p1d - p0d);

					//toa
					double rtoa = rref / mu0;
					double gtoa = gref / mu0;
					double btoa = bref / mu0;
					double ntoa = nref / mu0;

					//dem
					float demval = 0.0;
					{
						double currlon = leftLon + resolution * ix;
						double currlat = topLat - resolution * iy;
						if (currlon > 180.0) currlon = currlon - 360 ;
						int demx = (currlon - demTrans[0]) / demTrans[1];
						int demy = (currlat - demTrans[3]) / demTrans[5];
						if (demx >= 0 && demx < demdata.getXSize() && demy >= 0 && demy < demdata.getYSize())
						{
							demval = demdata.getValuei(demx, demy, 0);
							if (demval < 0) demval = 0;
						}
					}
					

					//atm result
					float sph[4], rho[4], tray[4], tog[4];

					//atm-corr
					int acok = getatmvariables(
						mu0, mu1, rph,
						demval,
						sph, rho, tray, tog
					);
					if (acok == 0)
					{
						//path correct
						if (u10tbb < 230)
						{
							for (int ib = 0; ib < 4; ++ib) rho[ib] *= 0.65;
						}
						else if (u10tbb < 290)
						{
							double tcorr = 0.65 + 0.35 * (u10tbb - 230) / 60;
							for (int ib = 0; ib < 4; ++ib) rho[ib] *= tcorr;
						}

						//6s sph
						sph[0] = 0.141; sph[1] = 0.10934; sph[2] = 0.0488;	sph[3] = 0.0157;
						tray[0] = 0.831 ; tray[1] = 0.8769; tray[2] = 0.947; tray[3] = 0.98;
						tog[0] = 0.99; tog[1] = 0.98; tog[2] = 0.944 ; tog[3] = 0.98;

						//ac
						double nac = correctedrefl(ntoa, tray[3], tog[3], rho[3], sph[3]);
						double rac = correctedrefl(rtoa, tray[2], tog[2], rho[2], sph[2]);
						double gac = correctedrefl(gtoa, tray[1], tog[1], rho[1], sph[1]);
						double bac = correctedrefl(btoa, tray[0], tog[0], rho[0], sph[0]);

						//limb vza blend
						double vzablend = 1.0;
						double nirfraction = 0.07;
						int t0t1 = MAX(t0, t1);
						if (t0t1 >= 8500)
						{
							vzablend = 0;
						}
						else if (t0t1 > 7500)
						{
							vzablend = 1.0 - (t0t1 - 7500) / 1000.0;
						}
						
						double nac2 = nac;// *vzablend;
						double rac2 = rac;// * vzablend;
						double gac2 = gac;// * vzablend;
						double bac2 = bac;// * vzablend;

						double newgac = (1 - nirfraction) * gac2 + nirfraction * nac2;
						
						//gamma
						int outr = pow(rac2 , gammainv) * 255;
						int outg = pow(newgac, gammainv) * 255;
						int outb = pow(bac2 , gammainv) * 255;

						outr = jac255(linear255(rac2 ) );
						outg = jac255(linear255(newgac)  );
						outb = jac255(linear255(bac2 )  );

						outr = MIN(MAX(0, outr), 255) * vzablend + (1 - vzablend) * nightr ;
						outg = MIN(MAX(0, outg), 255) * vzablend + (1 - vzablend) * nightg;
						outb = MIN(MAX(0, outb), 255) * vzablend + (1 - vzablend) * nightb;

						outdata.setValuei(ix, iy, 0, outr);
						outdata.setValuei(ix, iy, 1, outg);
						outdata.setValuei(ix, iy, 2, outb);
					}
					else
					{
 
						//bad ac
						outdata.setValuei(ix, iy, 0, nightr);
						outdata.setValuei(ix, iy, 1, nightg);
						outdata.setValuei(ix, iy, 2, nightb);
					}
				}
				else
				{
					//night
					outdata.setValuei(ix, iy, 0, nightr);
					outdata.setValuei(ix, iy, 1, nightg);
					outdata.setValuei(ix, iy, 2, nightb);
				}
			}
			else
			{//out space
				outdata.setValuei(ix, iy, 0, 0);
				outdata.setValuei(ix, iy, 1, 0);
				outdata.setValuei(ix, iy, 2, 0);
			}
		}

		int per1 = iy * 100 / ysize;
		if (per1 != per0)
		{
			per0 = per1;
			cout << per0 << ". "  ;
		}
	}

	cout << "100. start saving ... " << endl;
	double trans[6] = {leftLon, resolution, 0.0, topLat, 0.0, -resolution };
	outdata.copyTrans(trans);
	outdata.copyProj("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]");
	outdata.save(outfile);

	cout << "done." << endl;
	return 0;
}





